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We present an ab initio determination of the spin response of the unitary Fermi gas. Based on 
finite temperature quantum Monte Carlo calculations and the Kubo linear-response formalism, we 
determine the temperature dependence of the spin susceptibility and the spin conductivity. We show 
that both quantities exhibit suppression above the critical temperature of the superfluid-to-normal 
phase transition due to Cooper pairing. The spin diffusion transport coefficient does not display 
a minimum in the vicinity of the critical temperature and drops to very low values D s w 0.8 h/m 
in the superfluid phase. All these spin observables show a smooth and monotonic behavior with 
temperature when crossing the critical temperature T c , until the Fermi liquid regime is attained at 
the temperature T* , above which the pseudogap regime disappears. 

PACS numbers: 03.75.Ss, 05.60.Gg, 51.20.+d, 05.30.Fk 



There are basically two flavors of superfluids - 
fcrmionic and bosonic. The bosonic superfluid is real- 
ized when typically weakly interacting bosons condense 
and form a Bose- Einstein condensate (BEC). The typi- 
cal fermionic superfluid is of Bardeen-Cooper-Schrieffer 
(BCS) type, where fermions in time-reversed orbitals 
form weakly bound Cooper pairs and turn into a BEC of 
Cooper pairs. The BEC superfluidity vanishes when the 
condensate fraction ceases to have a macroscopic value 
with increasing temperature and the long-range coher- 
ence is lost. On the other hand, when a BCS super- 
fluid undergoes a phase transition to a normal state, 
Cooper pairs break due to thermal motion and there are 
no more bosonic constituents left to form a BEC. With 
the experimental confirmation of the BCS-BEC crossover 
in fermionic cold gases [H, 0], it became clear that one 
can have a new kind of system where both bosonic and 
fermionic superfluid features are present at the same 
time. The paradigmatic example is the unitary Fermi 
gas (UFG) where, unlike the BCS or BEC cases, the in- 
terparticle interaction is strong, and the binding energy 
of the Cooper pair is comparable to the Fermi energy. 
It is believed that the unpolarized unitary Fermi gas, in 
a temperature region just above the critical temperature 
T c , exists in a state which is neither fully bosonic nor 
fully fermionic in character, called the pseudogap state, 
widely studied in high-T c superconductors (HTSC) 
This is a temperature regime where a significant fraction 
of the Cooper pairs is present, even though the long- 
range coherence among them is lost, as is superfluidity. 
While the existence of the pseudogap regime in HTSC 
is an experimentally well known fact, the nature of the 
corresponding regime around the critical temperature T c 
of the UFG has remained an open question. It is a tanta- 
lizing question, whether the pseudogap regime exists in 
dilute neutron matter as well. The properties of the neu- 
tron superfluid in the neutron star crust are very similar 




FIG. 1: (Color online) Temperature evolution of the density 
of states profiles, extracted from the QMC simulations [3. 
The (blue) lines marked T* ~ 0.2 ef correspond to the onset 
of the pseudogap regime, the (red) line marked T c — 0.15 ef 
corresponds to the critical temperature. The profiles for T — 
0.21 ef and T = 0.12 ef are compared with the density of 
states for a Fermi liquid (dashed lines). 



to those of the unitary g ctS, clS highlighted by Bertsch's 
1999 Many-Body X challenge. The reliable calculation 
of the neutrino processes in neutron stars, controlled by 
the neutron spin response [Hj , is a long-standing problem 
and the present results will likely shed new light on these 
phenomena. 

The most striking feature of the pseudogap regime is 
the behavior of the fermionic density of states, which 



2 



shows a dramatic depletion at the Fermi level. This was 
confirmed to exist in quantum Monte Carlo (QMC) sim- 
ulations of the UFG jg, 0], see Fig. 1, as well as in ex- 
periments The transition from a pseudogap to a 
normal state in the UFG is somewhat similar to a gas- 
plasma transition, where no discontinuities are observed, 
which makes it difficult to observe. Indeed, no observable 
imprints on the thermodyna mic prop erties have been de- 
tected in experiments so far 111. Jl2j. while at the same 
time theory also predicts none [a, lD, |l3j, [l4| . A direct 
measurement of the local density of states of a UFG in 
a trap is desirable. However, one can suggest different 
kinds of measurements as well. A strongly paired, but 
not necessarily superfluid, system would respond qualita- 
tively differently to an external probe than a non-paired 
system, if one were to try to separate the two fermions 
in a pair. In particular, in a paired system the spin sus- 
ceptibility and spin conductivity should be significantly 
suppressed when compared to the unpaired regime. One 
should therefore observe a marked difference between the 
response of a system in the pseudogap regime from a nor- 
mal Fermi liquid one, which is the expected behavior at 
temperatures greater than T*. 

Inspired by recent measurements of various spin re- 
sponses like the spin susceptibility [Hf, spin trans por t 
coefficients [3] or dynamic (spin) structure factors [17j . 
we present here an ab initio evaluation of the spin suscep- 
tibility Xs and the spin conductivity a s for unpolarized 
homogeneous unitary Fermi gas. We show that both of 
these quantities carry a strong signal indicating survival 
of the Cooper pairs above the critical temperature. Our 
new results are consistent with previous studies @, 0], 
where the existence of the pseudogap regime was found 
between the critical temperature T c — 0.15(l)£i? and 
T* = 0.19(2) £p, where £p = p F /2m is the free Fermi 
gas energy, and pf = fi.(37r 2 n) 1 / 3 is the Fermi momen- 
tum corresponding to the total particle number density 
n. While in the previous estimation of the T* temper- 
ature finite resolution of the method provided only the 
lower bound, the analysis of the spin responses allows us 
to provide a significantly more accurate interval for the 
pseudogap onset. 

Moreover, the computed responses allow us to ex- 
tract the temperature dependence of the spin diffusion 
coefficient from first principle calculations. There has 
been considerable speculation as to whether the trans- 
port coefficients possess universal lower bounds imposed 
by quantum mechanics. The best known example is a 
conjecture formulated by Kovtun, Son, and Starinets of 
the existence of a lower bound r\j s ^ h/(Airk B ) on the 
ratio of the shear viscosity r\ to the entropy density s 
for all fluids [18j . While in our previous work [l9[ we 
found that path integral Monte Carlo (PIMC) calcula- 
tion is compatible with a well-defined minimum for the 
r)j s ratio in the vicinity of the critical temperature, here 
we show that the spin diffusion does not exhibit a similar 



minimum as a function of temperature. 

To determine the spin properties of the UFG we em- 
ploy the PIMC technique on the lattice, which provides 
numerically exact results, up to quantifiable systematic 
uncertainties (for details see Ref. [13J ) . Since the pseudo- 
gap regime is expected to exist in a rather small temper- 
ature region 0.15 < T/ef ^ 0.2 it needs to be checked if 
it survives when the thermodynamic V — » oo and contin- 
uum n — > limits are recovered. This step is of great im- 
portance as the critical temperature approaches the value 
T c = 0.15(1 ) Ef only in the thermodynamic and contin- 
uum limits [131 ] . At the same time the temperature T* 
above which there are no Cooper pairs left, and which re- 
flects short range correlations among particles, does not 
show a similar strong volume dependence. To check the 
stability of the results, as the thermodynamic and con- 
tinuum limit are approached, we performed simulations 
using three lattices N x — 8, 10, 12 with corresponding av- 
erage density n ^ 0.08, 0.04, and 0.03, respectively. The 
systematic errors related to finite volume effects as well 
as effective-range corrections, are estimated to be likely 
~ 10% - 15%, while the statistical errors of the PIMC 
data are below 1%; see [2(| for more details. Henceforth, 
we define units h = m = feg = 1- 

The spin susceptibility as well as the spin conductiv- 
ity can be theoretically determined using linear response 
theory via the Kubo relations. The uniform static spin 
susceptibility \ s = d(n^ — n,)/d((X^ — //.) is obtained 



as 
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where s z q = h q ^ — n q ^ represents a difference between 
spin-selective particle number operators in Fourier rep- 
resentation h qX = J2 p a{(p)a\(P + q), P = 1 / T 
is the inverse temperature and (. . .) stands for the 
grand-canonical ensemble average. The imaginary-time 
dependence of an operator is generated as 0(t) = 
e T(H-nN)Q e -r(H-nN)^ where fj ig the Hamiltonian of 

the system, \x is the chemical potential, and N is the 
particle number operator. The expectation values can 
be evaluated directly for q = 0. In this case the spin 
operator s q=0 commutes with the Hamiltonian and the 
expectation value is r independent. Consequently, the 
QMC calculation consists of the evaluation of the expec- 
tation value of a two-body operator, and the static spin 
susceptibility can be computed very accurately within 
our framework. 

In Fig. [SJ the static spin susceptibility x s m units of 
free Fermi gas susceptibility Xo = 3rt/2£p is shown for 
temperature range 0.1 < T/sp ^ 0.5. The results on 
8 3 , 10 3 , and 12 3 lattices exhibit satisfactory agreement 
with each other and no systematic trend in the data has 
been detected as we approach the thermodynamic and 
the continuum limit. For temperatures 0.25 Ef— 0.5 Ef no 
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FIG. 2: (Color online) The static spin susceptibility as a 
function of temperature for an 8 3 lattice solid (red) circles, 
10 3 lattice (blue) squares and 12 3 lattice (green) diamonds. 
The vertical black dotted line indicates the critical tempera- 
ture of the superfluid to normal phase transition T c = 0.15 ef- 
For comparison, the Fermi liquid theory prediction and recent 
results of the T-matrix theory produced by Enss and Hauss- 
mann [2^] are plotted with solid and dashed (brown) lines, 
respectively. The experimental data point from Ref. [ill is 
also shown. 



FIG. 3: (Color online) The spin drag rate Y s d = n/a s in 
units of Fermi energy as a function of temperature for an 8 3 
lattice solid (red) circles, 10 3 lattice (blue) squares and 12 3 
lattice (green) diamonds. Vertical black dotted line locates 
the critical temperature of superfluid to normal phase tran- 
sition. Results of the T-matrix theory are plotted by dashed 
(brown) line [2g]. The inset shows the extracted value of the 
contact density as function of the temperature. The (purple) 
asterisk shows the contact density from the QMC calculations 
of Ref. [32] at T = 0. 



strong temperature dependence of the spin susceptibility 
is observed. The latter is well below the susceptibility 
of the free Fermi gas with a value around x s ~ 0A5 % , 
which is in qualitative agreement with the Fermi liquid 
picture as well as results of other groups 22H26j . In the 
interval T* — 0.20 — 0.25 Ef we find the beginning of 
pronounced suppression of the spin susceptibility, which 
we associate with the existence of Cooper pairs in the 
system. Note that already at T c the spin susceptibil- 
ity is about half its value at the onset of suppression 
(roughly T — 0.25 £f)- Thus, two temperature scales are 
clearly distinguishable: the critical temperature of the 
superfluid-to-normal phase transition T c — 0.15 Ef, and 
the onset of the Cooper-pair formation T*. 

The static spin conductivity o~ s represents another 
quantity which is expected to be strongly affected by 
the presence of the Cooper pairs as it measures re- 
sponse of the spin current j s = j-f — j± once the weak 
external F force which couples with opposite signs to 
the two spin populations is applied to the system, i.e., 
j s = <j s F '. In order to extract the spin conductivity we 
consider the Kubo formula, which relates the frequency- 
dependent spin conductivity to the corresponding spec- 
tral density: cr s (ui) — nps\q — 0,u>)/uj; while the 
static spin conductivity is defined in the limit of zero 
frequency: a s = lim w _j. + o~ s (ui) . The spectral density 
Ps (?i w ) is related to the imaginary-time (Euclidean) 



current-current correlator provided by the PIMC method 
1 
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([f qt (r)-f qi (T)][f_ qt (0)-f_ qi (0)]), 



by inversion of the spectral relation 



(2) 



where j qX (r) stands for the third component of Fourier 
representation of the spin-selective current operator. To 
invert Eq. (j3]) we have applied the methodology which 
combines two complementary methods: singular value 
decomposition and maximum entropy method, both de- 
scribed in Ref. 27j. An additional a priori informa- 



tion includes the non-negativity of the spin conductivity 
o~ s (uj) ^ 0, a Lorentzian-like structure at low frequen- 
cies (Drude model), and the asymptotic tail behavior 
a B (uj — > 00) = C/(37rw 3 / 2 ), where C is the Tan contact 
density [25[ . The contact density was extracted from the 
tail of the momentum distribution n(p) ~ Cp~ 4 , using a 
technique similar to that of Ref. [28(. As the universal 
decay of the tail distribution starts around p/pf ~ 2, 
we were unable to extract the contact for the N x = 8 
lattice. In that case we used the value of C extracted 
from the N x = 10 lattice. The inset of Fig. [3] shows 
the temperature dependence of the contact density used 
as a priori information. For calculations in the N x = 12 
lattice, we found that the signal-to-noise ratio for the cor- 
relators at T < 0.16 ef is insufficient to perform a stable 
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FIG. 4: (Color online) The spin diffusion coefficient obtained 
by the Einstein relation D s = a a jxa as function of tempera- 
ture. The notation is identical to Fig. O 



reconstruction of the spectral density. For more details 
about the reconstruction process see the Supplemental 
Material 

In Fig. [3] we present the inverse of the static spin con- 
ductivity called the spin drag rate T s d = n/o~ s , which 
is the rate of the momentum transfer between fermions 
with opposite spins. For T = 0.25 Ef — 0.5 Ef, no strong 
temperature dependence is observed. For all three lat- 
tices the spin drag rate exhibits a significant enhance- 
ment above T c , in the interval T* = 0.20 — 0.25 ef, which 
is consistent with the occurrence of the spin susceptibil- 
ity suppression. Such an enhancement is expected for a 
system with strong correlations between particles of op- 
posite spins. 

Finally, the computed spin susceptibility and spin con- 
ductivity allow us to extract the spin diffusion coefficient 
D s in a fully ab initio manner. In the hydrodynamic 
regime it defines the proportionality between the spin 
current j s and spatially varying polarization by Fick's 
law j s = — _D s V(n-f — The spin diffusivity D s can 
be related to the spin conductivity and the spin sus- 
ceptibility by the Einstein relation D s — a s /xs- In 
Fig. |4] we show the temperature evolution of the spin 
diffusion coefficient. In the normal phase, for temper- 
atures 0.25 Ef — 0.5 Ef the diffusivity is approximately 
constant D s ss 1.8. Surprisingly, we find that the spin 
diffusion coefficient decreases substantially when the sys- 
tem enters the pseudogap regime, acquiring eventually a 
value around D s ss 0.8 in the superfluid phase. Such a 
low value can be understood as a quantum limit for this 
transport coefficient. The bound originates from kinetic 
theory, where D s ~ vl, v is the average particle speed, 
and I is the mean free path. For a strongly correlated 
system, the product of v ~ pf ~ n 1 / 3 and I ~ rt -1 / 3 
cancels the density dependence, giving D s ~ 1. 

In Ref. [29J the existence of a minimal value of the dif- 



fusivity was predicted for a temperature somewhat be- 
low the Fermi temperature, within Landau-Boltzmann 
theory. Recently, it was reported that Luttinger-Ward 
theory sets the minimum D s ^ 1.3 at T = 0.5 Ef [HI]. 
Our ab initio calculations do not confirm the presence 
of a minimum for the spin diffusion coefficient down to 
T = 0.1 Ef, and they do not rule out possibility that 
the diffusivity D s decreases further when temperature is 
lowered. 

Our results for the spin susceptibility and the spin drag 
rate deviate from the recent measurements of the MIT 
group [l6| extracted from fully polarized cloud collisions. 
This technique, in which two noninteracting clouds col- 
lide, and in which pairs do not exist and likely are not 
formed, is less suitable to probe the low temperature 
regime, where pairs already exist and their presence is 
of crucial importance. Pair formation in such an exper- 
iment would require three-body collisions. Recent the- 
oretical simulations [13, 0, [3l[ demonstrated that this 
experiment can be explained assuming that the measure- 
ment explores a noncquilibrium state associated with a 
quasirepulsive Fermi gas and two-body collisions alone. 
On the other hand, the technique based on speckle imag- 
ing of spin fluctuations, used to measure the spin sus- 
ceptibility of the system in thermal equilibrium (ljj . is 
in remarkable agreement with our theoretical results; see 

Fig. m 

In summary, we have presented results for the spin 
response of the UFG at finite temperature, obtained 
through an ab initio PIMC approach. The spin sus- 
ceptibility and the spin conductivity bear signatures of 
Copper-pair formation above the critical temperature 
T c - 0.15 £ F uptoT* w 0.20-0.25 e F . The spin diffusion 
coefficient does not display a minimum in the vicinity of 
the critical temperature, but instead drops to very low 
values D s s» 0.8 in the superfluid phase. We showed that 
the spin response of a unitary Fermi gas is not affected 
by the superfluid to normal transition, but only by the 
presence of Cooper pairs, and all these spin observables 
show a smooth and monotonic behavior up to the tem- 
perature T*, where the pseudogap cease to show up in 
the density of states. 
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Note added. Similar studies of Xs(T) and density of 
states in an attractive 2D Hubbard model (in the context 
of cuprates) have been performed in Refs. 



331. There 



are, however, qualitative differences between the physics 
of the attractive 2D Hubbard model and dilute Fermi 
gases in two dimensions, which are practically noninter- 
acting; see Chap. 7 in Ref. [2j. 
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SYSTEMATIC ERRORS 

The results of lattice simulations are affected by two 
types of errors: statistical and systematic. The former 
are explicitly under control as we collected ensemble of 
10 4 uncorrelated samples at each temperature, reducing 
the statistical uncertainties below 1%. The systematic 
errors arise due to the finite volume of the box V = N 3 
(we set lattice spacing to unity) and the finite density n 
of the system. The latter is related to corrections coming 
from the nonzero effective range r c g and exclusion of the 
universal high momenta tail in the occupation probability 
n(p) ~ Cp~ 4 . In the continuum limit n — > both correc- 
tions are negligible as pi?r e ff — > and p max /pF — > +oo, 
where p max = t is the maximum momentum available 
on the lattice. In our simulations we have reached val- 
ues PFr c tf — 0.54, 0.43, 0.39 and p max /p F i= 2.4, 3.0, 3.3, 
respectively for lattices N x = 8, 10 and 12. In par- 
ticular the moderate value of p m axlPF leads to signifi- 
cant bias of the current-current correlator . It is 
easy to note that contribution from the universal spin 
conductivity tail a s {to — > oo) = C / (Sttuj 3 ^ 2 ) leads to 
Gi n) {r = 0) = +oo. This condition is violated in PIMC 
simulations. In Fig. [1] we present the evolution of the 
Gs (t) correlator for temperature T — 0.42 Ep when the 
continuum limit is approached. The influence of the cor- 
relator by the systematic errors for low imaginary times 
is clearly visible. Note however, that the most valuable 
information about the low frequency part of the spin con- 
ductivity comes from the values of the correlator located 
around t//3 = 1/2, where the contribution from the uni- 
versal tail is marginal. 
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FIG. 1: (Color online) The dimensionless current-current 
correlator for temperature T — 0.42 ef as function of the 
imaginary time for 8 3 lattice solid (red) circles, 10 3 lattice 
(blue) squares and 12 3 lattice (green) diamonds. The inset 
shows violation of the spin conductivity sum rule in the lattice 
simulations as the function of temperature. Solid lines are a 
guide to the eye. 



In Refs [J Q it was shown that the spin /-sum rule 
implies (m = 1) 



+OC i 

— a s (uj) = ri. 



It is straightforward to prove that 

r+oo j,., r P 



f °°— cx s (u;) = f G^(r)dr, 



(1) 



(2) 



and thus the sum rule (JTJ) can be checked directly in 
the lattice simulations and used as an approximate mea- 
sure of the systematic errors. The deviation Sn = 
n — Jq Gs*^ (t) dr from the expected result is presented 
in the inset of the Fig. [TJ We have found that the dis- 
crepancy does not exceed 15% for each case, and the 
strongest violation of the sum rule is in the low temper- 
ature regime. This is in agreement with the previous 
estimations of the systematic errors reported in Ref. 
Moreover, we observe that the sum rule is typically un- 
derestimated (Sn > 0). This fact as well as the finite 
value of the correlator Gs (r = 0) are due to the ex- 
istence of the cut-off in the momentum space which im- 
plies the existence of the cut-off in the frequency space 
at uj max /e F « (p m a X /p_F) 2 . 

Our simulations predict the temperature of Cooper- 
pair formation to be T* /ep s» 0.20 — 0.25, which does 
not change significantly with the box size and the density 
(no systematic trend has been detected). On the other 
hand, the critical temperature T c does show finite-range 
effects. One may therefore ask the following question: 
Can the systematic errors merge the two distinct tem- 
peratures T c and T* into one temperature? Since the 



temperature of pair formation is not visibly influenced 
by the systematic corrections, we assume that T* is not 
lower than 0.2 ef in the continuum limit. In order to shift 
the critical temperature T c = 0.15 Ef to T* m 0.2 ef, 
one should assume that the systematic errors generated 
by effective-range corrections are at least 30%, which is 
far above our estimations. Moreover, results reported in 
Refs. (see also Ref. Q) do not reveal such strong 

departure of T c from its continuum limit value for den- 
sities considered in this work. A more precise analysis, 
utilizing improved actions and operators such as those 
constructed in Refs. d, 0] is required to rigorously quan- 
tify the uncertainties, but this is beyond the scope of this 
work. 

Recently, Castin and Werner [l(| made the claim that 
on a lattice, the manner the ultraviolet cutoff is im- 
plemented might lead to corrections which do not van- 
ish when the continuum limit is taken, especially if one 
adopts a so-called spherical momentum cutoff. 

In lattice calculations one introduces a single-particle 
propagator with an ultraviolet momentum cutoff 



Gi(r,t) 



00 duj_ 

-co 27ri 
2 1.2 , 



D 



d 3 k exp(ifc ■ r 
(2^)3 



iujt) 



+ i0+ - e(fc) 



(3) 



where e(fe) = H k /2m, r 



i y ,n z )l stand for dis- 



crete lattice points and where I is the lattice constant, the 
domain of integration D is either \k\ < A = tt/1 (spherical 
momentum cutoff) or \k x \ < A, \k y \ < A, \k z \ < A (cubic 
cutoff). If the lattice is infinite n XyVyZ run over all integer 
numbers. In the case of a finite box with periodic bound- 
ary conditions n x ,y,z run over a finite interval only and in 
this case the integral becomes a discrete sum in which d 3 k 
is replaced by (2ir) 3 /V, where V is the volume of the box. 
From the single-particle propagator one then constructs 
the two-particle propagator G\(ri — r2,t)G\(ri — T2,t) 
assuming that the two particles have to be at the same 
lattice point in order to interact. In the energy momen- 
tum representation this two-particle propagator reads 



G 2 (K,E) 
d 3 q 



1 



(2tt) 3 E + M+- e{K/2 + q)- e{K/2 - q) 



•(4) 



where D\2 = D\ PI Z?2 is the intersection of the two- 
momentum domains defined above either for a cubic or 
spherical cutoff for each single-particle momentum K /2± 
q respectively, and K is center of mass momentum of 
the pair. Usually the strength of the on-site interaction 
between two particles is obtained by assuming that both 
E and K vanish in the Lippmann-Schwinger equation, 
namely (we use convention V(r\ — r>i) = —god(ri — ra)) 



1 

9o 



d 3 q 



1 



iA 



ATrh 2 a J D (2tt) 3 2e(q) 4irh 2 a 2ir 2 h 2 



(5) 



Notice than in this case Dn = D. In the case of a cu- 
bic momentum cutoff the coefficient in front of the last 
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term is slightly different [10( . If the interacting pair has 
a finite center of mass momentum K and total energy 
E = h 2 K 2 /4m + H 2 k 2 /m, the corresponding Lippmann- 
Schwinger equation reads 



1 



1 



1 



g t{K, E) K 2 J Dl2 (2tt) 3 k 2 + i0+ - q 2 



m / 1 



4irh 2 \a 



ik 



mA ikm 



2ir 2 h 2 4irh 2 



in 



2n 2 h 2 



o 



Q±dq± I ,., 



dq x 



^h?- q \~\K\l2 k 2 + i0+ -q±-qi 



m mA m\K\ 
4nh 2 a ~ 2ir 2 h 2 + 8tt 2 H 2 



(6) 



where t(K, E) is i-matrix of the pair. We neglected terms 
of order k 2 and K 2 in evaluating the above integrals. Now 
we introduce an on-site coupling constant g for each value 
of the center of mass energy/momentum 



1 m mA 
g ~ Airh 2 a 2-K 2 h 2 
m mA 



Airh 2 a 2ir 2 h 2 



1 



m\K\ 
8ir 2 h 2 

\K\ 
4A 



(7) 



In a unitary Fermi gas the average value of | K | is of the 
order of kp , and thus this amounts to a correction of the 
order of kp/A, which vanishes in the continuum limit. 
Notice that by defining in this manner g the pole of the 
t- matrix is always at k = i/a irrespective of the center of 
mass momentum of the pair. At unitarity one has 



9 = 9o 1- 



\K\ 
4A 



5o 



90 n (hL 
1 A 



(8) 



thus a relatively small and systematic correction to the 
coupling constant, which can be safely evaluated in per- 
turbation theory and which vanishes in the continuum 
limit, when A — > oo. In our simulations we neglect this 
correction to the coupling constant and the results corre- 
spond to an interaction strengths shifted slightly towards 
the BCS limit. 

Another lattice artifact is the undesired finite range of 
the order of the lattice constant, which the contact inter- 
action acquires. There have been in literature a number 
of attempts to remedy this deficiency, one of the most 
popular being the modification of the dispersion relation 
of the kinetic energy in various ways, see Ref. and a 
number of references therein to earlier works. These sug- 
gestions amount to modifying the single-particle kinetic 
energy used in the two-particle propagator Q such that 



e(k) 



h 2 k 2 h 2 k 2 



2m 



2m 



k 2 

A 2 



(9) 



where the parameters a, A and possible higher order 
terms are chosen so as to either minimize the effective 



range, or even make it exactly zero, when such a dis- 
persion relation is used for the two-particle propagator 
(H]) in the Lippmann-Schwinger equation. In particular, 
for a > the wave function of the two-particle bound 
state would have the form [exp(— br) — exp(— cr)]/r (with 
b > and c > 0) instead of the expected exp(— r/a)/r, 
and the 1/r singularity at the origin disappears. The 
negative effect of such a procedure is that the expecta- 
tion value of the kinetic energy of the gas is modified 
as well and the properties of the free Fermi gas are not 
reproduced correctly in the corresponding QMC calcu- 
lation. Since the single-particle momentum distribution 
at large momenta does not satisfy anymore the expected 
behavior n(k) cx C/A; 4 all Tan relations 



111 are violated 



in the entire BCS-BEC crossover interval at all temper- 
atures. The energy of the interacting system would not 
satisfy anymore the expected high-temperature behavior 
E(T) 3TN/2 + . . . for T — y oo and the virial expansion 
would be violated as well. 



ANALYTIC CONTINUATION & INVERSION 

To compute the frequency dependent spin conductivity 
one has to invert the integral equation, for zero momen- 
tum transfer q = 0: 

/•OO 

G ( s J) (t)= pi" } (uj)K(uj,T)duj, (10) 
Jo 

with the kernel 

cosh^(r-/3/2)] 
K{U ' T)= sinhM/2) ' (U) 

Analogously to the extraction of the shear viscosity (see 
in our previous work [IH), it is convenient to rewrite the 
Eq. (fTU)) in the form: 



G 5 ® {n 



n JO 



<j s (uj) K{w, r) dui, 



(12) 



where K(lo,t) = luK(uj,t) and a s (uj) = npi^^ 
Note that we have introduced explicitly the frequency 
cut-off (see discussion in the previous section). The cor- 
relator G^i^ (t) is provided by the PIMC method for 
at least 51 points in imaginary time r, uniformly dis- 
tributed in the interval [0, f3}. As the correlator sat- 
isfy Gs (t) = G^J^ (f3 — t) we restrict the inversion to 
the interval r 6 [0,/3/2], and apply the symetrization: 
G [ s jj \r) <h- [G^'V) + G { r\p - r)]/2 in order to de- 
crease statistical fluctuations. 

The inversion is based on combination of two comple- 
mentary methods: Singular Value Decomposition (SVD) 
and self-consistent Maximum Entropy Method (MEM), 
described in details in paper (l~3| . The inversion proce- 
dure consists of two steps: 



8 




FIG. 2: (Color online) The frequency dependent spin con- 
ductivity obtained from the Quantum Monte Carlo calcula- 
tions with 10 3 lattice for three different temperatures: upper 
panel T = 0.15 sf (at the critical temperature), middle panel 
T = 0.20 £f (pseudogap regime) and lower panel T = 0.30 ef 
(normal phase). The gray bands show uncertainty of the re- 
construction. In the insets by (blue) points the corresponding 
PIMC correlators are presented. The line shows quality of re- 
production of the PIMC data. 



1. the SVD solution a s (uS) is created. This step ex- 
ploits only the information contained in the corre- 
lator G { J j) {r). 

2. the self-consistent MEM produces final solution 
o~ s (uj), which is forced to satisfy external con- 
straints of the form: 



P[aM(w)]=*M(u>), 



,0'i) 



c 



3nw 3 / 2 ' 



for lu > ui taih 



(13) 
(14) 



where P stands for the projection operator onto the 
SVD subspace. To represent the MEM solution, 
we used a mesh of 120 points uniformly distributed 
within interval [0,u; ma x]. 

The self-consistent MEM needs to be supplemented by a 
class of a priori models for the solution. The class has 
been defined as 



M(w,{j[i,7,c,a x ,a2}) = f(u>, {a x , a 2 }) 



C 



37TW 3 / 2 



+ [!-/(<*>, "2})] if*, 7, c}), (15) 



where 



and 



/(w,{o! x ,a 2 }) = e Q i c 



1 



- 1 



1 + e Q i( CJ_Q 2) ' 



7T (iO 



Ai) 2 +7 2 ' 



(16) 



(17) 



The parameters {/i, 7, c, a x , a 2 } describe admissible de- 
grees of freedom of the model and are adjusted auto- 
matically in the self-consistent manner. To initialize the 
model for the first iteration we fit it to the SVD solution. 

The algorithm includes three "bootstrap" parameters: 
i) Wmax - the upper limit of the integration, ii) w tail - 
the point where the universal tail behavior is imposed, 
iii) a - the parameter which governs the relative impor- 
tance of the a priori model to the data. The bootstrap 
sample consists of about 200 launches of the algorithm 
with randomly generated (from some reasonable inter- 
vals) bootstrap parameters. Since each solution satisfy 
the integral equation (fT2|) also the solution produced as 
an average over the collection satisfies it as well. This 
solution is regarded to be the most favorable scenario, 
while the standard deviation over the bootstrap sample 
is associated with the uncertainty of the inversion proce- 
dure. 

We have used the Lorentzian-likc structure C(oj) at low 
frequencies (Drude model) for the conductivity as it cor- 
rectly reproduces T-matrix results for u) < Sf as reported 
in [H . Moreover we have checked the stability of the solu- 
tion when one uses the Gaussian-like structure instead of 
the Lorentzian. We have not found significant differences 
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for the static spin conductivity a s as we change the low 
frequency part of the model. 

In Fig. [5] we present a sample of results for three 
temperatures obtained for 10 3 lattice. In the normal 
phase T > 0.25 ep we observe that a s {uj) possesses a 
Lorentzian-like structure at low frequencies with a max- 
imum in the vicinity of zero frequency. Once we enter 
the pseudogap regime, the suppression of <J s (uS) devel- 
ops for the zero frequency, due to the onset of pairing 
correlations. 
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